//  time_series.cpp
//
//      Copyright 2009 Stefan Fruhner <stefan.fruhner@gmail.com>
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
// MA 02110-1301, USA.


//#define DEBUG

#include <iostream>
#include <cstdlib>
#include <vector>
#include <set>

// includes from mylibs
#include <mylibs/mymesh.hpp>
#include <mylibs/mymatrix.hpp>
#include <mylibs/mystring.hpp>
#include <mylibs/mymesh.hpp>
#include <mylibs/lookuptable.hpp>
#include <mylibs/myinifiles.hpp>
#include <mylibs/cmdline.hpp>
#include <mylibs/gipl.h>
#include <mylibs/myio.hpp>

#include <myImg/IGBheader.hpp>

using namespace std;
using namespace mylibs;
using namespace mylibs::cmdline;

#ifndef print
//#define print(x) (cout << __FILE__ << ":" << __FUNCTION__ <<"() (" <<__LINE__ << ") : " << x << endl)
#define print(x) (cout << " - " << x << endl)
#endif

#ifndef printe
#define printe(x) (cerr << __FILE__ << ":" << __FUNCTION__ <<"() (" <<__LINE__ << ") "<< "Error : " << x << endl)
#endif

// assign some value a to a variable b, does conversion and print the new
// result to stdout
#ifndef ASSIGN
#define ASSIGN(a,b) { a >> b; TELL(b); }
#endif

class Voxel{
	vector<unsigned int> index; // holds the x,y and z-index

	public:
		Voxel(unsigned int x, unsigned int y=0, unsigned int z=0){
			index.push_back(x);
			index.push_back(y);
			index.push_back(z);
		}

		unsigned int& operator[](size_t i){
			return this->index[i];
		}

		void info(){
			if (index.size() == 0 ) return;
			cout <<"[ ";
			cout << index[0] << ", ";
			cout << index[1] << ", ";
			cout << index[2] << " ]" << endl;
		}
};

bool verbose = false;
Point pt1(0.,0.,0.);
Point pt2(0.,0.,0.);
bool pt1_set = false;
bool pt2_set = false;
vector<unsigned int> voxel1;
vector<unsigned int> voxel2;
set<size_t> skip_t;
bool continuation_mode 	= false;
vector<mystring> igbfiles;
vector<mystring> meshes;
double tol = 1.e-99;

mystring xlabel = mystring("");
mystring ylabel = mystring("");
mystring title  = mystring("");

void write_gnuplot_file(mystring datfile){

	if (not IO::is_executable("gnuplot")){
		cmdline::warning("gnuplot was not found. I'd prepare some graphs if you install it.");
		return;
	}

	mystring gnufile(datfile.file_base()     + ".gnu");

	cout << " - Writing gnuplot file " << gnufile << endl;
	mystring epsfile(datfile.file_base(true) + ".eps");
	mystring fn2(datfile.file_base(true) + ".dat");
	ofstream gnu;
	gnu.open(gnufile.c_str());


	gnu << "set terminal postscript eps enhanced 20 dl 4. lw 2.0\n"
		<< "# set pointsize 2.0\n"
		<< "# set terminal x11\n" << endl
		<< "set output \""<< epsfile <<"\"\n" << endl
		<< "set xlabel \""<< xlabel <<"\"\n" << endl
		<< "set ylabel \""<< ylabel <<"\"\n" << endl
		<< "#set xrange [0.:500.] \n"
		<< "#set yrange [-1.:1.] \n"
		<< "set format y \"%1.4f\"\n";

	if (pt2_set or voxel2.size() > 0)
		gnu	<< "plot \""<< fn2 << "\" u 2:($4-$3) ";
	else
		gnu	<< "plot \""<< fn2 << "\" u 2:3 ";


	if (title.empty()){
		gnu << "notitle w l" << endl
			<< "# pause -1" << endl;
	} else {
		title.replace("_", "\\_");
		gnu << "title '" << title << "' w l" << endl
			<< "# pause -1" << endl;
	}

	gnu.close();

	// run gnuplot
	string cmd("gnuplot " + gnufile.file_strip_path());
	cmdline::run_command(cmd.c_str(), datfile.file_path().c_str());
	return;
}

template <typename T>
bool is_zero(T *data, size_t items){
	for (size_t i = 0; i < items; ++i){
		if (data[i] != 0.) return false;
	}
	return true;
}

template <typename T>
int read_IGBfile(IGBheader& h, ostream &o, int item){

	int items = h.x()*h.y()*h.z();
	T *data = new T[items];
	int slicesize = items * h.data_size();

	int numread = slicesize;

	if (not continuation_mode){

		do{
			int i_t = (gztell( (gzFile) h.fileptr()) - 1024)/ slicesize;

			double t = h.org_t() + i_t * h.inc_t();
			numread = gzread( (gzFile)  h.fileptr(), data, slicesize); 	// read one slice
			if (numread != slicesize) break; 					// if no data was read
			if (is_zero(data, items) ) {
				cmdline::warning("Skipping value at t = " + toString(t) + ". All values are zero !");
				continue;
			}

			// check if this slice shall be skipped
			if (skip_t.find((size_t)i_t) != skip_t.end()) continue;

			size_t idx1 = h.index(voxel1[0], voxel1[1], voxel1[2]);		// get the 1D index of the voxel in the mesh
			o << i_t <<"\t"<< t << "\t" << data[idx1];					// write it to the stream
			if (voxel2.size() > 0){
				size_t idx2 = h.index(voxel2[0], voxel2[1], voxel2[2]);	// get the 1D index of the voxel in the mesh
				 o << "\t" << data[idx2];								// write it to the stream
			}
			o << endl;
		}while (numread == slicesize);

	}
	else { // continuation mode

		SurfaceMesh mesh(meshes[item], verbose);
		if (item == 0 and not pt1_set) { // if a voxel-ID was given get the coordinates
			pt1 = mesh.p[voxel1[0]]; // we need the coordinates of the point

			do{
				int i_t = (gztell( (gzFile) h.fileptr()) - 1024)/ slicesize;
				double t = h.org_t() + i_t * h.inc_t();
				numread = gzread( (gzFile) h.fileptr(), data, slicesize);
				if (numread != slicesize) break; 								// if no data was read
				if (is_zero(data, items) ) {
					cerr << "Warning: Skipping value at t = " << t << ". All values are zero !" << endl;
					continue;
				}
				// check if this slice shall be skipped
				if (skip_t.find((size_t)i_t) != skip_t.end()) continue;

				T &value = data[voxel1[0]];						// get the value
				o << i_t <<"\t"<< t << "\t" << value << endl;	// write it to the stream
			} while (numread == slicesize);

		}else { // if the user specified coordinates directly
			vector<double> D;
//			size_t interpolation_tet  = mesh.find_element(pt, D, tol);
			size_t interpolation_tet  = mesh.find_closest_element(pt1, D, tol);
//			mesh.fprintf_barycentric_coordinates(stdout, D);
			if (interpolation_tet >= mesh.elements() )
				cmdline::exit("Could not find surrounding tet. Error occured at mesh "
								+ toString (item) +". You may want to increase -tol");

			do {
				int i_t = (gztell( (gzFile) h.fileptr()) - 1024)/ slicesize;
				double t = h.org_t() + i_t * h.inc_t();
				numread = gzread( (gzFile) h.fileptr(), data, slicesize); 		// reading complete slice
				if (numread != slicesize) break; 								// if no data was read
				if (is_zero(data, items) ) {
					cerr << "Warning: Skipping value at t = " << t << ". All values are zero !" << endl;
					continue;
				}

				// check if this slice shall be skipped
				if (skip_t.find((size_t)i_t) != skip_t.end()) continue;
				// barycentric interpolation of the function value within this tet
				T valueBC = (T) 0.;
				for (size_t i = 0; i < D.size()-1; i++){
					valueBC += D[i+1] * data[mesh.f[interpolation_tet].v[i]];
				}

				// simple mean value in tet
//				T value = (T) 0;
//				for (size_t i = 0; i < mesh.f[interpolation_tet].size(); i++){
//					value += data[mesh.f[interpolation_tet].v[i]];
//				}
//				value /= mesh.f[interpolation_tet].size();

				o << i_t <<"\t"<< t << "\t" << valueBC  << endl;	// write data to the stream
			} while (numread == slicesize);
		}
	}
	delete[] data;

	if (verbose){cmdline::ok();}
	return numread;
}
int main(int argc, char** argv){
	mystring inifilename("");

//!####################### Load settings from file #############################"

	myIniFiles ini(argc, argv);
	if (not inifilename.empty())
		print(" - Reading settings from ini file " << inifilename);

	ini.set_info("time_series ... Extrahieren von Zeitreihen aus CARP igb.gz\
 - files und Verwerfen von Zwischenschritten, um files zu verkleinern.");

	ini.register_param("meshes", "i", "mesh file(s) in CARP or Tetgen\
 style", true);
	ini.register_param("out", "o", "Folder where the output is written to hard disk.");

	ini.register_flag("verbose", "v", "Additional information on stdout.");

	ini.register_param("igbfiles", "d" , "filename(s) of the igb files\
 to load", true);

	ini.register_param("voxel", "V" ,"<int int int> voxel to compute\
	the time series for. The 3 ints define one voxel in the mesh given by\
 x,y and z index. For irregular meshes - where 3 indices make no sense -\
 mostly the x-index is used, whereas all others are 0.");

 ini.register_param("voxel2", "U" ,"<int int int> second voxel that is\
  used for difference plots (NOT implemented in continuation mode).");

 ini.register_param("point", "p" ,"<double, double double> coordinate in the\
 mesh, that will be tracked in continuation mode. In normal mode the closest\
 voxel is selected", true);

 ini.register_param("point2", "2" ,"<double, double double> coordinate in the\
 mesh.  In normal mode the closest voxel is selected and the difference signal\
 f(p2)-f(p1) is computed.", true);

	ini.register_flag("continuation_mode", "c", "continuate the physical\
 coordinates of the voxel if more than one mesh is given. If not set \
 in all meshes the value at the given voxel ID is used");

	ini.register_flag("gather", "g", "Data from multiple igbfiles is\
 gathered in one output file.");

	ini.register_param("skip_t", "s", "[List of integers] Times to skip from the anaylsis,e.g. ");

	ini.register_param("xlabel", "x", "x label of the gnuplot graph.");
	ini.register_param("ylabel", "y", "y label of the gnuplot graph.");
	ini.register_param("title" , "T", "Title of the gnuplot graph."  );
	ini.register_param("tol",    "t", "tolerance value , usally a small number. Default: 1.E-5");
	ini.check(3);

	verbose = ini.exists("verbose");
	tol 	= ini.read("tol"   , tol   );
	xlabel  = ini.read("xlabel", xlabel);
	ylabel  = ini.read("ylabel", ylabel);
	title 	= ini.read("title" , title );
	mystring output = ini.read("out", mystring());

	ini.readVector("meshes"	,	meshes	 );
	ini.readVector("igbfiles",	igbfiles );

	ini.readVector("voxel" ,		voxel1);
	ini.readVector("voxel2",		voxel2);

	vector<size_t> skip;
	mystring test = ini.read("skip_t", mystring());
	if (test.contains("-") ) {
		vector<mystring> t = test.splitv("-");
		size_t min = 0;	t[0] >> min;
		size_t max = 0;	t[1] >> max; max++;
		TELL(min); TELL(max);
		for (size_t i = min; i < max; i++)skip.push_back(i);

	}	else{
		ini.readVector("skip_t", 		skip);

	}
	skip_t = std::set<size_t>(skip.begin(), skip.end());
//	info(skip);

	continuation_mode = ini.exists("continuation_mode");
	bool gather_data  = ini.exists("gather");

	vector<double> point;
	ini.readVector("point", point);

	vector<double> point2;
	ini.readVector("point2", point2);

	// grep a given point
	if (point.size()) {
		if (point.size() != 3) cmdline::exit("Points must have 3 coordinates !");

		pt1.set(point[0],point[1],point[2]);
		pt1_set = true;
		pt1.info();

		if (voxel1.size() > 0 ) {
			info(voxel1);
			cmdline::exit("You have to decide if you search for point or voxel. Both is not allowed");
		}


	} else if ( voxel1.size() > 0 ){ // if a voxel was given
		// fill indices to 3 items
		for (size_t i = voxel1.size(); i < 3; ++i) voxel1.push_back(0);
	}

	if (point2.size()) {
		if (point2.size() != 3) cmdline::exit("Points must have 3 coordinates !");

		pt2.set(point2[0],point2[1],point2[2]);
		pt2_set = true;
		pt2.info();

		if (voxel2.size() > 0 ) {
			info(voxel2);
			cmdline::exit("You have to decide if you search for point or voxel. Both is not allowed");
		}
	} else if ( voxel2.size() > 0 ){ // if a voxel was given
		// fill indices to 3 items
		for (size_t i = voxel2.size(); i < 3; ++i) voxel2.push_back(0);
	}

//! ########################## Error handling ##################################
	if (voxel1.size() > 3 or voxel1.size() > 3 )
		exit("Voxels must have 3 indices for x,y and z-direction!");

	if ( igbfiles.size() == 0 ) cmdline::exit("Need igb files to run ... " );
	else
	if (meshes.size() > 0 and meshes.size() != igbfiles.size()){
		cerr << "Need " << meshes.size() << " igbfiles, but "
			 << igbfiles.size() << "given." << endl;
		exit(-1);
	}

	if (ini.exists("point") and meshes.size()== 0)
		exit("If you define a point, then you will have to define a mesh");

	if (continuation_mode and meshes.size()== 0){
		continuation_mode = false;
		cmdline::warning("Continuation mode can only be used when \
you define a mesh. Resetting continuation_mode=false.");
	}

	if (not gather_data and continuation_mode){
		cmdline::msg("Flag 'gather' was enabled due to continuation mode.");
		gather_data = true;
	}

//! ########################## create output path   ############################

	IO::mkdir(output);

//! ############################### main program ###############################

	// convert point into voxel coordinates of the closest voxel
	if (pt1_set or pt2_set){
		SurfaceMesh mesh(meshes[0]);
		if (pt1_set){
			voxel1.clear();
			voxel1.push_back(mesh.find_closest_node(pt1));
			voxel1.push_back(0);
			voxel1.push_back(0);
		}

		// convert point into voxel coordinates of the closest voxel
		if (pt2_set){
			voxel2.clear();
			voxel2.push_back(mesh.find_closest_node(pt2));
			voxel2.push_back(0);
			voxel2.push_back(0);

			if (continuation_mode) cmdline::warning("Attention : Point2 will\
 be ignored in continuation mode !");
		}
	}

	filebuf fb;
	string fn;
	for (size_t i = 0; i < igbfiles.size(); i++){ // gehe alle igb files durch
		if (verbose)
			cmdline::section(igbfiles[i].c_str());

		gzFile file = gzopen(igbfiles[i].c_str(), "r");
		if (not file) {
			cmdline::warning("File " + igbfiles[i] + " could not be opened.");
			if (i != 0)continue; else exit(-1);
		}

		IGBheader h(file);
		// simply compute the index to trigger an exception if needed
		h.index(voxel1[0], voxel1[1], voxel1[2]);

		gzseek(file,1024, 0); // seek to the end of the header
		// decide if a new file has to be opened
		if ( 		(not gather_data)  		 // open new file for each input
				or  (gather_data and i == 0) // gather data: open only one file
			) {
			mystring file(igbfiles[i]);

			// set the filename for output
			if (file.lower().endswith("igb.gz"))
				fn = igbfiles[i](0, -7).file_base(true);
			else
				fn = igbfiles[i].file_base(true);

			// append
			if (not output.empty())
				fn = output.path_join(fn);


			fn += "_"   + toString(voxel1[0]).fill(6,'0') + ".dat";
//					 	+ toString(voxel1[1]) + "_"
//					 	+ toString(voxel1[2]) + ".dat";

			cmdline::msg(" - Opening file " + fn);
			fb.open(fn.c_str(),ios::out);
		}

		if (verbose){
			TELL(h.x());
			TELL(h.y());
			TELL(h.z());
			TELL(h.t());
			print("x = "<< h.org_x() <<" + i_x * " << h.inc_x());
			print("y = "<< h.org_y() <<" + i_y * " << h.inc_y());
			print("z = "<< h.org_z() <<" + i_z * " << h.inc_z());
			print("t = "<< h.org_t() <<" + i_t * " << h.inc_t());
			TELL(h.dim_x());
			TELL(h.dim_y());
			TELL(h.dim_z());
			TELL(h.dim_t());
			TELL(h.unites_x());
			TELL(h.unites_y());
			TELL(h.unites_z());
			TELL(h.unites_t());
			TELL(h.fac_x());
			TELL(h.fac_y());
			TELL(h.fac_z());
			TELL(h.fac_t());

			info(voxel1," - Voxel 1");
			info(voxel2," - Voxel 2");
		}

		ostream os(&fb);
		switch (h.type()){
			case IGB_BYTE:  read_IGBfile<unsigned char >(h, os, i); break;/* -- byte ------------*/
			case IGB_CHAR:  read_IGBfile<char  >		(h, os, i); break;/* -- Char ------------*/
			case IGB_SHORT: read_IGBfile<short >		(h, os, i); break;/* -- short -----------*/
			case IGB_LONG:  read_IGBfile<long  >		(h, os, i); break;/* -- long ------------*/
			case IGB_FLOAT: read_IGBfile<float >		(h, os, i); break;/* -- float -----------*/
			case IGB_DOUBLE:read_IGBfile<double>		(h, os, i); break;/* -- Double ----------*/
			case IGB_INT:	read_IGBfile<int   >		(h, os, i); break;/* -- integer ---------*/
			case IGB_UINT:	read_IGBfile<unsigned int  >(h, os, i); break;/* -- unsigned integer-*/
			case IGB_COMPLEX:	/* -- 2 x float (real part, imaginary part) - */
			case IGB_D_COMPLEX:	/* -- 2 x Double (real part, imaginary part)  */
			case IGB_RGBA:		/* -- 4 x byte (red, green, blue, alpha) ---- */
			case IGB_STRUCTURE:	/* -- Structure ----------------------------- */
			case IGB_POINTER:	/* -- void * -------------------------------- */
			case IGB_LIST:		/* -- List   -------------------------------- */
				printe("Datatype not implemented" << h.type()<< " " << h.TypeID2Name(h.type()));
				break;
			default:
				printe("Datatype unknown : " << h.type() << " " << h.TypeID2Name(h.type()));
				exit(-2);
				break;
		}

		// determine if the file needs to be closed
		if (not gather_data or i == igbfiles.size()-1 ){
			fb.close();

//			mystring xl(h.unites_t());
//			if (xl == "mv") xl = "mV";
//			xl = "time / " + xl;
			if (xlabel.empty())
				xlabel = "time / " + mystring(h.unites_t());
			if (ylabel.empty())
				ylabel = "data / " + mystring(h.unites() );
//			mystring yl(h.unites()); yl = "{/Symbol F}_e / " + yl;
			write_gnuplot_file(fn);
		}

		gzclose(file);
	}

//! ############################### finalization ###############################
	print("Saving settings ...");
	ini.writeVector("meshes",	meshes	);
	ini.writeVector("igbfiles",	igbfiles);

	if (continuation_mode)
		ini.write("continuation_mode", continuation_mode);

	if (pt1_set)
		ini.writeVector("point",	point);
	else
		ini.writeVector("voxel",	voxel1);

	if (pt2_set)
		ini.writeVector("point2",	point2);
	else
		ini.writeVector("voxel2",	voxel2);

	ini.save();

	return 0;
}
